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ABSTRACT 

The density of stars in galactic bulges is often observed to be flat or slowly rising inside the influence 
radius of the superniassive black hole (SMBH). Attributing the dynamical friction force to stars moving 
more slowly than the test body, as is commonly done, is likely to be a poor approximation in such 
a core since there are no stars moving more slowly than the local circular velocity. We have tested 
this prediction using large-scale A^-body experiments. The rate of orbital decay never drops precisely 
to zero, because stars moving faster than the test body also contribute to the frictional force. When 
the contribution from the fast-moving stars is included in the expression for the dynamical friction 
force, and the changes induced by the massive body on the stellar distribution are taken into account, 
Chandrasekhar's theory is found to reproduce the rate of orbital decay remarkably well. However, this 
rate is still substantially smaller than the rate predicted by Chandrasekhar's formula in its most widely- 
used forms, implying longer time scales for inspiral. Motivated by recent observations that suggest 
a parsec-scale core around the Galactic center SMBH, we investigate the evolution of a population 
of stellar-mass black holes (BHs) as they spiral in to the center of the Galaxy. After ~ 10 Gyr, we 
find that the density of BHs can remain substantially less than the density in stars at all radii; we 
conclude that it would be unjustified to assume that the spatial distribution of BHs at the Galactic 
center is well described by steady-state models. One consequence is that rates of capture of BHs by the 
SMBH at the Galactic center (EMRIs) may be much lower than in standard models. When capture 
occurs, inspiraling BHs often reach the gravitational-radiation-dominated regime while on orbits that 
are still highly eccentric; even after the semi-major axis has decreased to values small enough for 
detection by space-based interferometers, eccentricities can be large enough that the efficient analysis 
of gravitational wave signals would require the use of eccentric templates. We finally study the orbital 
decay of satellite galaxies into the central region of giant ellipticals and discuss the formation of 
multiple nuclei and multiplet of black holes in such systems. 

Subject headings: black hole physics-Galaxy:center-Galaxy:kinematics and dynamics-stellar dynamics 



1. INTRODUCTION 

Dynamical friction plays a central role in many as- 
trophysical contexts. It drives t he orbital inspiral 
and merger of satellit e galaxies (e.g. Murai & Fuiimoto 
119801 : llbata fc Lewisi [l998: van den Bosch ct al. 199a) 
and the form a tion of massive black hole b i naries 
(e.g. lOuinla nI 119961: iMilosavlievic fc Merrittl 120011 : 
iMakino & Fu nato I l20M r" and it is the fundamental 
mechanism leading to mass segregation in dense stellar 
systems (e.g. [Bahcall fc Wolf 1977; Frcitag et al. 2006; 
iHopman fc Alexander II2OO60 . 

Chandrasekhar formulated the principle of dynami- 
cal friction under the assumptions of an infinite, ho- 
mogeneous and isotropic field of stars (Chandrasekhar 
1943). Despite these simplifications, his theory has 
been shown to work remarkably well in a wide va- 
riety of more general situations. Dynamical friction 
can be understood as the drag induced on a test par- 
ticle by the overdensity (i.e., the gravitational wake) 
that is raised behind it by the deflection of star s 
(|Danbv fc Camml[l957l : IKalnais I [l97l IMulder I Il98l . 
The surprisingly good agreement between theory and 
numerical results may be attributed to the fact that 
the wake is a local structure, and over small spatial 
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scales, the stellar background appears nearly homoge- 
neous ([Weinberg I [l986.) . On the other hand, numeri- 
cal studies have revealed a few, astrophysically impor- 
tant contexts in which Chandrasekhar's theory appears 
to break down. These include the deceleration of a ro- 
tating stellar bar ([Weinberg 1985), inspiral in harmonic 
(constant-density) cores (Hernandez fc Gilmore 199^; 
IGoerdt et an[2006l: iRead et al. |[2006l: llnoue 1 12009) , and 
the o rbital evolution of a displac ed supermassive black 
hole (jGualandris fc Merritt |[2008D . 

In this paper, we present a comprehensive study of 
dynamical friction in the nuclei of galaxies containing a 
dominant central point mass. In particular, we investi- 
gate the case of shallow density profiles around supermas- 
sive black holes (SMBHs). Such nuclei appear to be com- 
mon, and perhaps even generic. For instance, the lumi- 
nosity profiles of br ight elliptical galaxies a lways exhibit 
flat c entral cores (Fcrrarcs e et al. I 119941 : iLauer et al. I 
119951 ). Even the Milky Way, which was long believed 
to have a steeply-rising mass density near Sgr A*, is 
now believed to have a p arsec-scale core (iBuchholz et al. I 
I2OO9I : IDo et al. I [2OO9I : IBartko et al.1 mM . Similar 
models may also be applicable to dark matter halos, 
if the centr al point mass is identifled with the stel- 
lar spheroid (iBorriel lo fc Sal ucci |[200lHBinnev fc Evans I 
[200I iSpekkens et al. , 2005, ). 

Theoretical treatments of dynamical friction make a 
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surprising prediction about the frictional force in such 
systems. Essentiahy ah of the decelerating force is pre- 
dicted to come from stars that are moving more sfowly 
than the test body. But the phase-space density of a 
galaxy with a shallow density cusp around a SMBH falls 
to zero at low energies: below a certain radius (roughly 
1/2 the core radius), there are no stars locally that move 
more slowly than the circular velocity at that radius. 
Chandrasekhar's formula, in its most widely-used form, 
would predict no frictional force. In the case of an eccen- 
tric orbit that passes in and out of the core, the frictional 
force would be small near periapsis, leading to a rapid in- 
crease in orbital eccentricity - the opposite of the usual 
assumption. Our numerical experiments reveal that the 
frictional force does not drop precisely to zero in such 
nuclei. We show that the evolution can be well described 
by a more general form formula that includes a contribu- 
tion to the force from stars moving faster than the test 
mass. In this sense, our results affirm the correctness of 
Chandrasekhar's physical picture, but only if the proper 
field-star velocity distribution is used (as opposed to, say, 
a Maxwellian), and only if the usual simplifying assump- 
tions that lead to a neglect of the contribution of the fast 
stars to the frictional force are relaxed. 

In §2 we review Chandrasekhar's derivation of the dy- 
namical friction force and highlight the approximations 
that lead to the neglect of the contribution from the 
fast-moving stars. We also briefly discuss alternative 
treatments of dynamical friction. In §3 we use Chan- 
drasekhar's formulae to integrate the equations of motion 
of a massive body and follow its inspiral into the center 
of a model designed to represent the Galactic center. In 
§4 we use large-scale A^-body simulations to test the the- 
ory in the case of inspiral of massive objects in a nuclear 
star cluster with a flat density profile. §5 investigates the 
formation of the gravitational wake in the self consistent 
simulations. Applications of our results to a variety of 
astrophysical problems are discussed in §6, and §7 sums 
up. 

2. DYNAMICAL FRICTION 

The motivation for the A^-body experiments described 
in this paper is the existence of physically interesting 
models of galactic nuclei in which the standard dynam- 
ical friction formula predicts little, or zero, frictional 
force. We begin in this section by re-deriving the stan- 
dard formula, noting the simplifying approximations that 
are usually made. We then present the more general form 
of Chandrasekhar's formula that includes contributions 
from field stars of all velocities, not just those that move 
more slowly than the test body at infinity, and we evalu- 
ate the expected contribution from the fast-moving stars 
in our models. We also compute how the fast- and slow- 
moving stars contribute differently to the steady-state 
density wake, using a technique first applied by Mulder 
(1983). Finally we comment on perturbative approaches 
to computing dynamical friction that relax the assump- 
tion of an infinite homogeneous medium. The results ob- 
tained in this section constitute a set of baselines against 
which the A^-body results can be compared. 

2.1. Chandrasekhar's treatment 

iCha ndrasekhar ' (1943') derived the coefficient of dy- 
namical friction by summing the encounters of a test 



body with passing stars, assuming that the unperturbed 
motion of the test body was linear and unaccelerated, 
and that the field star distribution was infinite and ho- 
mogeneous spatially and isotropic in velocity space. 

The velocity change of a test body of mass M in one 
encounter with a field star of mass m <C M is 
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where V is the relative velocity at infinity, p is the impact 
parameter, and po = GM/V"^ . The velocity change in 
equation ([Ij is parallel to the initial, relative velocity V 
before the encounter. In order to derive the coefficient 
of dynamical friction, one sums the velocity changes in 
the direction of motion of the test mass, per unit interval 
of time, over all impact parameters and over all values 
for the relative velocity at infinity. The summation over 
impact parameters, at fixed V , is achieved by multiplying 
equation ([1} by 2'KpnVdp, with n the number density of 
field stars, and integrating dp: 
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Under the assumption that A 
be written as 
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Terms beyond the first in brackets, the so-called "non- 
dominant" terms, are usually neglected. 

Returning to the more general form ([2]), the dynam- 
ical friction coefficient is obtained by a second integra- 
tion over field star velocities v^,. The relative velocity is 
V = V — Vi,, with V the velocity of the test star. Since 
equation 1^ gives the velocity change in the direction of 
the initial relative motion, it must be multiplied by 
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to convert it into a velocity change in the direction of 
the test star's motion, assumed here to be along the x 
axis. Letf{v^,)dv^, be the number density of field stars in 
velocity increment Vi,,Vi, + dVi,, normalized to unit total 
number. The dynamical friction coefficient is 
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where p = mn. 

Henceforth we assume that the field star distribution 
is isotropic in velocity space. Following Chandrasckhar] 
(1943), we represent the velocity-space volume element 
in terms of Vi, and V using 
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The result is 
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(The quantity J defined in equation (26) of Chan- 
drasekhar (1943) is equal to Sv-^T-l.) The integral that 
defines H turns out to have an analytic solution; the ex- 
pression is complicated and we do not reproduce it here. 
Chandrasekhar (1943) gave several approximate forms 
for H valid for Pmax/po ^ 1, e.g. his equation (30): 
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In the standard approximation (e.g. lRosenbluth et al. I 
[1957), the non-dominant terms are set to zero, and the 
velocity dependence of the logarithmic term in the inte- 
grand of equation ^ is ignored. Instead, one writes 
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and the lower bound Pmin is set to GM/v1^„ 
weighting function % then takes on the simple form 
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and the coefficient of dynamical friction is 



(Avii) ^ -AttG'^ M p X Air dv^ (^'^Y f{v^). (10) 

Equation ([TU]) reproduces the well-known result that only 
field stars with Vi, < v contribute to the frictional force. 

In this paper, we consider models for galactic nuclei in 
which the number of stars moving more slowly than the 
test body can be vanishingly small. In such models, one 
expects that a significant fraction of the frictional force 
might come from stars with Vi, > v. 

The distribution of field-star velocities in our models 
has the following form within the core: 



where the normalizing constant 
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corresponds to unit total number. This expression is 
equivalent to equation ([S]); it gives the local distribu- 
tion of velocities at a radius where the circular veloc- 
ity is Vc = {GM,/r)^^^ , assuming the density of field 



stars follows r The phase space density is zero for 

V* > Vcsc = 2^/'^Vc- 

Of more interest here is the behavior of / at small 
values of w*, and when 7 < 3/2; for such values of 7 
the phase space density diverges at = 2^Pvc- As 
7 1/2, the velocity distribution becomes progressively 
narrower, and in the limit, f{Vi,) is a delta-function at 
Vi, = 2^Pvc', in other words, all stars have zero energy. 
This may be seen as a consequence of the well-known 
fact that p oc r"*^-^ is the shallowest power law density 
profile consistent with an isotropic velocity distribution 
in a point-mass potential. 

In the case of a test body moving in a circular orbit 
with V = Vc, the number of field stars with v^, < v will 
drop as 7 approaches 1/2, and will equal zero in the 
limiting case 7 = 1/2. The standard dynamical friction 
coefficient, equation (|10l) . predicts zero frictional force in 
this limit. 

In this situation, it is clearly of interest to compute 
the contribution of the fast-moving stars to the total fric- 
tional force. We did this by evaluating H in its "exact" 
form, equation (pc)) . Figure [T] shows the results. In ad- 
dition to 7, the results depend on the parameter 

lnA = ln ( 1 (13) 



GM, 

which plays the role of Coulomb logarithm. We note the 
following results. 

• For 7 > 3/2, the contribution to the frictional force 
from the fast-moving stars is negligible, particu- 
larly when In A is also large. 

• For 7 < 3/2, the fast-moving stars contribute a 
progressively larger fraction of the total frictional 
force, particularly when In A is small. 

• When 7 — 0.55, near the limiting value, the total 
frictional force is small, and almost all of it comes 
from stars with v^, > v. 

• Whereas the contribution to the force from the 
slow-moving stars depends strongly on 7, the con- 
tribution from the fast-moving stars is almost in- 
dependent of 7. 

According to equation (O, the contribution of the fast 
stars must tend to zero as In A is made sufficiently large. 
This is consistent with Figure [TJ however, for 7 « 0.5, 
the value of In A required for the slow stars to dominate 
is far greater than any physically reasonable value. 

2.2. Mulder's treatment 

The foregoing treatment highlighted the contribution 
of the fast-moving stars, v^, > v, to the total frictional 
force. However it did not provide much insight into why 
the two populations contribute in such a different way 
to the force. Of course, the iV-body experiments de- 
scribed in this paper include both populations of stars. In 
the simulations, the field stars quickly establish a nearly 
steady-state distribution in a frame m oving with the test 
mass - a "dynamica l -friction wake" (|Danbv fc Camm~l 
[T95I IKalnais I [1971 : IMulder I [l98l . The over-density 
in the wake is responsible for the decelerating force that 
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Fig. 1. — Contribution to the total dynami cal friction force from stars moving faster, or more slowly, at infinity than the test body, 
assuming the velocity distribution of equation lllll l. The test body is assumed to be moving at the local circular velocity Vc- In these plots, 
the configuration-space density p remains fixed as 7 is varied. 



acts on the test body. A large fraction of the mass in 
the wake must be contributed by the fast stars, particu- 
larly in the case that the fast stars dominate the density 
at large distances. Why then do these stars contribute 
relatively little to the frictional force? 

One way to address this question is via the tech- 
nique of iMuldeFI ()1983[ ). Mulder computed the steady- 
state distributions of stars around a moving test mass, 
making e ssentia l ly the same assumptions as made by 
fChandrasekhar I (|1943| ). He did this by invoking Jeans 's 
theorem in a frame moving with the test mass, and show- 
ing that an isotropic f{Vi) at infinity could be expressed 
in terms of two of the integrals of motion in the Ke- 
pler problem. This then allowed him to compute the 
steady-state density, in the moving frame, at all loca- 
tions around the test mass. The dynamical friction force 
followed from a second integration of the density over 
space; Mulder showed that the results for the frictional 
force so obtained were consistent with Chandrasekhar's 
predictions, if Pmax were associated with the maximum 
dimension of the spatial grid used to carry out the force 
integration. 

Mulder's technique can be modified, to compute the 
separate contributions to the dynamical friction wake of 
the fast {vi, > v) and slow {vi, < v) stars; here, as above, 
Vi, refers to the field-star velocity at infinity. The results 
are shown in Figure [2l for 7 = 5/4. For this choice of 
7, the fast stars dominate the total density at infinity. 
The density that they generate near the test body is 
also higher, everywhere along the symmetry axis, than 
the density due to the slow stars. However the shapes 
of the two density wakes are very different: in the case 
of the fast stars, the wake is elongated counter to the 
direction of the test body's motion, while in the direction 
parallel to the motion, the change in density between the 
upstream and downstream sides of the test mass is much 
less than in the case of the wake produced by the slow 



stars. These two differences are responsible for the small 
contribution of the fast stars to the total frictional force 
(Figure 1), in spite of the higher density of those stars at 
infinity and in the wake. 

Comparison of the upstream and downstream densi- 
ties in Figure 2 also suggests why the relative contribu- 
tion of the fast stars to the frictional force drops off with 
increasing In A in Chandrasekhar's treatment (Figure 1). 
At large distances from the test body, the wake produced 
by the fast stars is nearly symmetric; the greatest asym- 
metry is in the region near the test mass. The wake 
generated by the slow stars, on the other hand, main- 
tains its asymmetry much farther from the test body. 
Roughly speaking, the density far from the origin in Fig- 
ure 2 is produced by stars with large impact parameters, 
and so increasing pmax in Chandrasekhar's treatment cor- 
responds to more heavily weighting the contribution from 
the slow-moving stars. 

2.3. Perturbative treatments; inhomogeneous systems 

In treatments like Chandrasekhar's and Mulder's, the 
unperturbed trajectories consist of straight lines. In re- 
ality, both test and field stars follow non-rectilinear or- 
bits about the center of the galaxy. Chandrasekhar's 
theory might be expected to give approximately cor- 
rect results even in this case, as long as Pmax ^ Pmin, 
since over many decades in scale the orbits of the 
field stars will appear nearly rectilinear as seen by the 
test body. But given certain assumptions, perturba- 
tion theory can be used to more correctly compute 
the response of the orbits in a galaxy to the presence 
of a perturbing potential (iLynden-Bell fc KalnaisllT972l : 
iTremaine fc Weinberg II1984I : iRauch fc Tremaine ]ll99i. 
One finds that the net torque on the test mass is due 
to orbits near resonance, i.e. orbits for which the fre- 
quencies associated with the radial and angular motions 
satisfy a relation liLOr + h^e — h^t = where the li are 




z z z 

Fig. 2. — Dynamical friction wakes, computed via Mulder's (1983) technique, assuming equation JSj with 7 = 5/4 for the velocity 
distribution at infinity; the test mass is located at the origin and is assumed to be moving at constant velocity v = Vc, as in Figure 1. The 
top panels show contours of the density, in a plane that contains the test body's velocity vector; the left panel shows the total density, the 
middle panel shows the density contributed by the stars with ti* < at infinity, and the right panel shows the contribution from stars with 
Vi, > V at infinity. Black (solid) curves show the total response from the indicated stars; blue (dashed) curves show the part of the response 
that is symmetric with respect to z; red (dotted) curves show the asymmetric part (only on one side), which is responsible for the frictional 
force. The contours are spaced logarithmically in density and the contour spacing is different in the three panels. The lower panels show 
the density along the symmetry axis, i.e. along a line through the test body in the direction of its motion. Units are G = M = v = 1. 



integers and fit is the frequency of rotation of the test 
mass (assumed to be on a circular orbit). The accel- 
eration induced by the resonant orbits depends on how 
quickly the orbit of the test mass is evolving; if orbital 
decay is very slow, the influence of a single resonance can 
build up, invalidating the perturbative assumption, while 
if it is too fast, the assumption of near-stationarity is vio- 
lated. Furthermore, in a real galaxy (or A'^-body system) 
the frequency spectrum of the perturbing potential is not 
made up of sharp lines, but rather is broadened by the 
time dependence of the decaying orbit and by the finite 
age of the galaxy. 

Due to the computational complexity involved, appli- 
cations of this approach have so far been limited to bod- 
ies following circular orbits in simple (Plummer, scale- 
free) galaxy models, and the results have mostly been 
interpreted as corrections to the predictions of M ulder 
and Chandrasekhar. For instance. iWeinbergl (|1986l ) em- 
phasized the similarity in the structure of the wake as 
computed via the perturbation formulae and via Mul- 
der's approach. The main element that the perturba- 
tive treatments add is a quantitative estimate of the 
Coulomb logarithm. Not surprisingly, none of these stud- 
ies has attempted to relate the frictional force separately 
to the "fast" and "slow" stars as they appear in Chan- 
drasekhar's treatment; doing so would be an ill-defined 
problem since all stars are included, self-consistently, in 



the perturbative treatments. Nevertheless, as far as we 
can tell, comparisons with Chandrasekhar's theory are 
always made via equation (jlO[) . which ignores the fast- 
moving stars. 

A potentially important application of the perturba- 
tive methods is to cases where the assumption of local- 
ity is violated. For instance, a satellite that orbits just 
outside a galaxy, where the local density is zero, would 
experience no frictional force if the local properties of the 
background were assumed to hold everywhere; in reality 
it feels a force due to polarization of the orbits inside the 
galaxy iPalm er fc" Papaloizou 1985). The models con- 
sidered in this paper constitute a second case where the 
assumption of locality may be inappropriate, since some 
of the frictional force acting on a test mass orbiting in 
the core will come from stars outside the core, where 
f{v) has a different functional form, including (for in- 
stance) some slow-moving stars. In lieu of such a cal- 
culation (and in view of the difficulties associated with 
interpreting the results; e.g. [Weinberg ()2004[ )). an A^- 
body treatment seems a logical first step. As we will see, 
Chandrasekhar's formula, in its more general form, turns 
out to reproduce the A^-body results quite well. 



3. ORBITAL EVOLUTION BASED ON 
CHANDRASEKHAR'S FORMULAE 
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We are interested in the orbital evolution of a massive 
body as it spirals in toward the center of a galaxy that 
contains a supermassive black hole (SMBH). In subse- 
quent sections, we present results from large-scale, direct- 
summation A^-body simulations. As a basis for compar- 
ison, we present in this section the predictions of Chan- 
drasekhar's approximate formula. We represent the stars 
via a smooth, fixed potential and integrate the equations 
of motion of the massive body in the fixed analytic poten- 
tial including a term that represents the non-conservative 
contribution of dynamical friction. 

We base our model for the stellar density on the ob- 
served distribution of old stars at the Galactic center 
(GC) . Number coui its (iBuchholz et al. 1 120091: iDo et al. I 
[20091 IBartko et al. l[2010ir are consistent with a density 
that follows a broken power-law: 
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where a is a parameter that defines the transition 
strength between inner and outer power law s and is 
the scale radius. Following iMerritt I (|2010[ ). we adopt 
r-Q — 0.3pc , a — 4 and 7e = 1.8 as fiducial values. The 
central slope 7 was left as a free parameter. The normal- 
izing factor Po was chosen in such a way that for each 
value of 7, the corresponding density profile reproduces 
the coreless density model: 
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outside the core. This choice of normalizing constant 
gives a mass density at Ipc similar to what various au- 
thors have inferred fe.g. lOh et al. l[2009l ) and implies a to- 
tal mass in stars within this radius of ^ 1.6x lO^Mopc"^. 

Assuming equal-mass stars of mass m and an isotropic 
velocity distr ibution, the lo cal two-body relaxation time 
is defined as (jSpitzer II1987I ) : 



0.33cr3 
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where InA is the Coulomb logarithm and a is the isotropic 
velocity dispersion; the latter can be computed from 
Jeans's equation, 



p{r)a{ry = G / dr'r'-" [M, + M^{< r')] p{r'). (17) 

J r 

Here M, is the mass of the central SMBH that we take to 
be 4 X IO^Mq (IGhez et al. 1 120081 : IGiUessen et an[2009h 

and Mi,{< r) is the total mass in stars within r. The 
total stellar mass contained within the SMBH infiuence 
radius {r^h ~ 2.5pc) is M^,{< rbh) ~ lO^M©; assuming 
solar-mass stars, the two-body relaxation time at rbh is 
tr(rbh) « 2 X IQiOyr. 

3.1. Circular Orbits 

The frictional acce leration on a point par ticle of mass 
M and velocity v is (IChandrasekhar |[T943l) 
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where F{< v,r) is the fraction of stars at r that are 
moving more slowly than v. This is the standard expres- 
sion, derived by ignoring the velocity dependence of InA 
when integrating over the field-star velocity distribution 
and setting the non-dominant terms to zero. As a result 
of these approximations, the frictional force is produced 
only by field stars with velocities less than v. Although 
equation (|18p was derived under the assumptions of an 
infinite and homogeneous background of stars, it has 
been shown to work reasonably well even for more gen- 
eral stellar distributions (White 1983; Lin & Tre mainel 
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Just et al. \mm . 

For a massive particle initially located at rbh on a cir- 
cular orbit, the inspiral time in the power-law density 
profile of equation with 7 — 1.8 (i.e., the coreless 

model) is 



ifr sa 6 X lO'^yr 



2.5pc j (' 



2.5pcy VlOOkms-i 
1 X IO^MqX f 7 \ 



(19) 



independent of the mass of the field stars if ill >> m. 

Figure [3] plots the relaxation time as a function of ra- 
dius for the same model, assuming InA — 15, m — Mq 
and adopting different values for the inner density slope 
7. It turns out that the isotropic distribution function 
corresponding to the adopted density law becomes 
negative at certain energies for 7 < 0.6. For this reason, 
we consider in the following only models with 7 > 0.6. 

FigureOalso shows the evolution of a 2 x lO'^M© black 
hole on a circular orbit starting from a galactocentric 
distance of 2.5 pc and using In A = 7. The orbit was nu- 
merically integrated by solving the system of first-order 
differential equations 



r = v, V = -V(j) + 



(20) 



with (j){r) the total gravitational potential produced by 
the stars and the SMBH: 



(r) 



GM. 



Mr) 



GM, 



-AttG 



dr'r'^ 



p(r') 



dr'r'p{i') 



(21) 



The numerical integration was performed using a 7/8 
order Runge-Ku tta algorithm with a variable time-step 
(|Fehlberg | |1968f ) in order to keep the relative error per 
step in energy, in the absence of dynamical friction, less 
than a specified value (10~*). When dynamical fric- 
tion was included, we checked the integration accuracy 
through the quantity E + E^f with E the energy per 
unit mass and i?df the work done by dynamical friction 
along the trajectory. The accuracy in this case was of 
the same order of that found in integrations without dy- 
namical friction. Th e function F(< f , r) was evaluated 
using the expression (jSzell et al. II2005I) : 




1 

F{< v,r) = l-- / 
pJa 

v/V2 _i f v/V2 



,1 dp 



— tan 



, (22) 




r (pc) Time (yr) 

Fig. 3. — Left panel: relaxation time versus radius for models based on the density law of equation 1)1411 . Right panel: Orbital decay 
of a 2 X 10^A/q massive body starting from a radius of 2.5pc. Here we used InA = 7. In both panels, various values of the inner density 
slope 7 were considered: (0.6, 0.8, 1, 1.25, 1.5, 1.8). 
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Fig. 4. — Fraction of stars F(< fcirci'') moving more slowly 
than the local circular velocity as a function of radius for 7 = 
(0.6,0.8,1,1.25,1.5,1.8). When 7 = 0.6, F is close to zero for 
r Si O.lpc . Hence, the frictional force acting on a massive particle 
which moves on a circular orbit drops essentially to zero at this 
radius. 

where E = + </)(r). 

At all radii, the relaxation time is much longer than the 
time required for the massive particle to reach the core. 
What happens next depends on 7: the orbital decay can 
essentially stall when 7 is small (i.e., '-^ 0.6), or continue 
rapidly if 7 is larger. 

The explanation of this behavior can be found in Fig- 
ure|3]which plots the fraction of stars moving more slowly 
than the local circular velocity Vcirc as a function of ra- 
dius, for various values of 7. When 7 = 0.6, F(< Wcirc, t) 
approaches zero at Tgt ~ O.lpc and consequently the dy- 
namical friction force drops drastically at this radius (see 
equation [TB]). The stalling observed in the orbital evo- 
lution for this value of 7 is therefore a consequence of 
the lack of slowly-moving stars in the core. However, 
the inspiral always continues into the very center since 



^"(< Wcirc, f ) > everywhere. 

For 7 > 0.6, the time required for dynamical fric- 
tion to bring a Wi^Mq black hole into the center, start- 
ing from a galactocentric distance of a few parsecs, is 
shorter than the two-body relaxation time evaluated at 
the SMBH influence radius tr(''bh)- On the other hand, 
the dynamical friction force decreases with the mass of 
the inspiraling object, and for M < IQ'^Mq the in- 
fall timescale can sign ificantly exceed a Hubble time. 
iMerritt fc Szell I (|2006l) found that t,{rbh) is also approxi- 
mately the timescale over which gravitational encounters 
change an initial density profile into the Bahcall-Wolf 
form, i.e., p cx r^^ ''^. We conclude that for a black 
hole of mass M > 10"^ M0, inspiral will occur in a mass 
profile that is almost independent of time. However, for 
7 ~ 0.6, the time required to reach a distance ~ O.Olpc, 
is still comparable with the local relaxation time. This 
will result in a substantial evolution of the stellar back- 
ground during the orbital decay. 

3.2. Eccentric Orbits 

In the case of an isotropic distribution function f{E) 
describing a power law density profile around a SMBH, 
if the gravitational potential produced by the stars is 
ignored (i.e., E « —GM,/r\,\y), then 

' 8 V ttS r(7 - 1/2) m (^GM,f Wo / 

, , (23) 

with 00 = GM./rbh (IMerritt 1 120 ll) . For 7 < 0.5, /(E) 
is undefined and so 7 « 0.5 is the shallowest density 
profile consistent with an isotropic velocity distribution 
around a SMBH. In the case 7 — 1.5, equation ([25)1 shows 
that the distribution function is a constant {f{E) = /q). 
If one writes 

If" 4 
p{r)F{< v,r) = p{r) x -r^47r / dv^vlfo = -tt /qw^ 

(24) 
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Fig. 5. — Left panel shows the time dependence of the orbital eccentricity of a M = 2 X 10^ Mq black hole. In the right panel the 
orbital evolution is shown in the eccentricity-semimajor axis plane. The inner cusp slopes are 7 = (0.6, 1, 1.25, 1.5, 1.8). Initial apoapsis 
and periapsis distances were 2.5 and 0.35pc respectively and initial semi-major axis was a = 1.4pc. The integrations terminated either 
when the semi-major axis of the black hole was O.Olpc or at 10*yr for 7 = 0.6. 



it can be immediately seen that the product p{r)F{< 
v,r) in equation will be a function of v only (e.g., 
Just et al. 2011). Under these circumstances, the coef- 
ficient of dynamical friction will have only a weak de- 
pendence on radius through the Coulomb logarithm. It 
can be shown that, in this case, the eccentricity of a 
massive body will remain unchanged during its motion, 
while dynamical friction will either circularize the or- 
bit for 7 > 1.5 or make it more eccen tric for 7 < 1.5 
(|Quinlan Ill996t IGould fc Guillen |[2Ciol . 

To evaluate the eccentricity evolution of a massive par- 
ticle in response to Chandrasekhar's dynamical friction 
formula, a numerical treatment is necessary. We there- 
fore carried out numerical integrations of the set of dif- 
ferential equations ((20|) as described above, adopting as 
before equations p4)) and ((2T|) for the (fixed) stellar po- 
tential. 

Figure [5] shows the results for A/ = 2 x 10^ Mq. The 
massive particle was initially placed at r = 2.5pc with a 
tangential velocity of ~ 0.36z;circ- With this initial con- 
figuration the body penetrates the inner core after few 
obits. Different values of the internal slope 7, ranging 
from 1.8 to 0.6 were adopted. As a proxy for the instan- 
taneous orbital elements, we computed over each radial 
period the largest and the smallest distance from the ori- 
gin (i.e. the SMBH) and defined these as the apoapsis 
Tap and periapsis rpor respectively. The eccentricity and 
semi-major axis were then computed using the Keplerian 
expressions 



' ap 



' per 



' ap 



1 



(25) 



The figure reveals a complex behavior of eccentricity on 
time. For 7 < 1.5 we distinguish three regimes. In phase 
I, the eccentricity decreases (even for 7 > 1.5). The 
duration of this phase is shorter for shallower profiles. 
After reaching a minimum, the eccentricity then increase 
rapidly with time (phase II). Finally, in phase III, the 
eccentricity either continues to increase, but more slowly 
than in phase II, or remains constant for 7 = 1.5. 



This evolution can be understood by considering the 
changes of rap and rper with time. In phase I, the black 
hole periapsis is close to the core radius where the dif- 
ference between the density models is small. As a conse- 
quence, the eccentricity evolution is nearly independent 
of 7 and the orbits circularize. In phase II, rpor is well in- 
side the core where the smaller dynamical friction results 
in a rapid eccentricity increase. Finally, in phase III, the 
orbit lies entirely inside the core. As a consequence of 
the declining dynamical friction at rap the eccentricity 
growth slows down. As predicted, for 7 = 1.5, the eccen- 
tricity remains unchanged in this phase. 

These results show that, in the presence of a fiat 
(7 < 1) density profile, a second black hole found initially 
on an eccentric orbit can acquire very large eccentricities 
(< 1) before entering the regime where relativistic ef- 
fects become important. In Section 16.21 we discuss in 
more detail how very large eccentricities may modify the 
expectations for the GW signal from massive black hole 
binaries for proposed space-based interferometers. 

In the first phase, when the periapsis is still outside 
the core, the orbit evolves completely in the outer cusp 
(7e = 1.8). Evolution in this regime could lead to a 
rapid circularization before the black hole reaches the 
inner core. To quantify the amount of circularization 
in this phase we computed a further orbit in the model 
with 7 = 0.8, adopting initially a larger semi-major axis 
(a ^ lOpc) and a smaller eccentricity (e = 0.3). The 
results of this integration (Figure (5]) show that the ec- 
centricity reaches a minimum value, e « 0.15, and then 
increases rapidly reaching e « 0.3 at rpor = O.lpc. At the 
end of the integration the orbit retains therefore a sub- 
stantial eccentricity (~ 0.4), even though it was almost 
circularized at the beginning of phase II. 

4. TV-BODY SIMULATIONS 

The numerical integrations of equation (jl8p presented 
above predict that a massive body that spirals in to the 
center of a galaxy containing a SMBH, and a nuclear 
star cluster with flat (7 < 0.6) density profile, will stall, 
at a radius that is roughly the core radius. Moreover, 
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Time(yr) Radius (pc) 

Fig. 6. — Left panel : eccentricity evolution for a 2 X IO^Mq black hole in a model with 7 = 0.8. The initial apoapsis and periapsis of 
the orbit are 12 and 7pc respectively which give a semimajor axis a 9pc. Right panel: eccentricity versus semimajor-axis (black line), 
apoapsis (green line) and periapsis (blue line). 



its eccentricity is expected to increase steeply once the 
orbital periapsis lies inside the core. Here we use A'^-body 
simulations to test these predictions. 

4.1. Initial Conditions and Numerical Method 

In order to generate equilibrium A'^-body models of the 
GC region that extend self-consistently to the Sgr A* 
influence radius (rbh ~ 2.5pc) we used the truncated 
mass model 



p{r) 



r ' 









r 


1+ 











(7-7c)/" 



with truncation function 



sech(x) + cosh(a;) 



Cir/n), (26) 
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Fig. 7. — Density profiles of equation II26I I with 7 = 
(0.6, 1, 1.25, 1.5, 1.8), ro = 0.3pc, a = 4 and truncation radius 
rt = 1.2pc (vertical dashed line). The dash-dotted line gives the 
coreless model of equation II15I I. 



With this choice, the density falls off exponentially at 
large radii (i.e., r > rt) while for r <^ rt, where 
({x) w l — x'^/8, the model reproduces almost exactly the 
density of equation p^ . As above, we chose tq = 0.3pc, 
a = 4, 7e = 1.8 and po = 1-3 x lO'^M©. Monte-Carlo ini- 
tial positions and velocities were then generated by nu- 
merically solving equation (j22[) : we stress that the equi- 
librium models so produced include self-consistently the 
effects of the gravitational force from the stars. Figure [7] 
shows the truncated density profiles for different values 
of 7 and rt — 1.2pc. 

The initial conditions wer e evolved usiiig the direct- 
summation code (/iiGRAPE (Harfst e t al. II2007D which 
uses a fourth-order Hcrmite integrator with a predictor- 
corrector scheme and hierarchical time steps. The per- 
formance and accuracy of the code depend both on the 
time-step parameter rj and on the smoothing length e. In 
what follows, we set rj = 0.01 and e = 5 x 10~''pc. With 
these choices, energy conservation was typically of order 
0.1% over the entire length of the integration. Most of 
the A'^-body integrations were carried out on the 32-node 
GRAPE cluster at the Rochester Institute of Technology. 
In addition, a few were carried out in serial mode using 
a Tesla C870 graphics processing unit with SAPPORO, 
a CUDA library that emulates double-precisio n force cal- 
culat ions on single precision hardware ( Gaburov et. al."l 
l2009h . 

Table 1 gives the parameters of the A'^-body models. 
The initial distance of the secondary black hole is given 
by Tin while its initial orbital eccentricity is ein- The 
quantity r* is the radius at which the initial mass in 
stars equals M, the mass of the second black hole. All 
of our A^-body models had nn < rt, so that the orbital 
evolution is expected to be very similar to that in the 
corresponding non-truncated models. In order to study 
the dependence of the results on the secondary black hole 
mass we run simulations with a range of masses, M = 
(2000, 5000, 10000, 5OOOO)M0. Two cases with nonzero 
initial eccentricities (runs Gl and G2, with ei — 0.54) 
were also considered. 
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M 
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Al 


0.6 


230 


1.2 


5 


22 





1 


0.07 


6.7 


A2 


0.6 


130 


1.2 


5 


38 





1 


0.07 


6.6 


Bl 


0.8 


230 


1.2 


5 


22 





1 


0.06 


6.9 


B2 


0.8 


130 


1.2 


5 


38 





1 


0.06 


6.9 


C 


0.6 


80 


0.6 


5 


26 





0.5 


0.07 


6.3 


D 


0.6 


130 


1.2 


2 


38 





0.3 


0.05 




E 


0.6 


130 


1.2 


10 


38 





1 


0.10 


6.4 


F 


0.6 


130 


1.2 


50 


38 





1 


0.18 


4.8 


Gl 


0.6 


200 


1.2 


5 


25 


0.54 


1 


0.07 


6.9 


G2 


0.6 


100 


1.2 


5 


50 


0.54 


1 


0.07 


6.9 



4.2. The Coulomb Logarithm 

In Table 1 we report the values of the Coulomb log- 
arithm extracted from each iV-body integration. The 
value of InA was obtained by minimizing the quantity: 

n 

^[r,(^)-r'(^,lnA)]^ (28) 

i=l 

outside a galactocentric radius r > 0.3pc. Here, n is the 
number of A'^-body data points, (i) is the position of the 
black hole in the A^-body simulation at time t, and r'{t) 
is its position at the same time evaluated by means of the 
Chandrasekhar's formula ([TS]). Since analytical expres- 
sions are not available for the trajectory of an inspiraling 
black hole, in order to obtain the expected position r'{t) 
at any given time, we first solved numerically the equa- 
tions of motion I^U^ and then built a spline interpolant 
from the results of the integration. This procedure was 
applied only in the part of the orbit outside the core, 
where equation (|18p is able to describe accurately the 
black hole orbit. In this way, unlike in most previous 
studies, we could obtain an estimate of the Coulomb log- 
arithm without making any assumptions about the ve- 
locity distribution of the field stars (e.g., that it followed 
a Maxwellian distribution). 

Our simulations do not show any obvious dependence 
of InA on either the number of particles or on the ini- 
tial eccentricity. We found an average value of InA = 
6.5 ± 0.2, in essentially perfect agre ement with the value 
reported bv lSpinnato et al. I (|2003[ ): InA ^ 6.6 ± 0.6. 

4.3. Results 
4.3.1. Circular Orbits 

The first simulations we performed consisted in evolv- 
ing the massive body on a circular orbit with initial ra- 
dius O.lpc (i.e., smaller than the stalling radius when 
7 ^ 0-6) and for a time corresponding approximately to 
300 orbits (i.e., ~ 4 x lO^yr at this distance). We used 
A^ = 130,000, M = 5OOOM0 and 7 = (0.6, 1, 1.5, 1.8). 
We also implemented a high-resolution simulation with 
N = 500,000 for the model with 7 = 0.6. As in most 
of the longer simulations of Table 1, the truncation ra- 
dius was rt — 1.2pc. These shorter integrations allowed 
us to study dynamical friction, while limiting the devia- 
tions of the models from their initial configuration that 
was found to occur on longer timescales as a result of 



two-body relaxation and perturbations from the massive 
object (see below) . The eccentricity of the orbit remained 
small during these integrations (e < 0.1). 

Figure [8] shows the time evolution of the semi-major 
axis of the orbits and the rate of orbital decay s — ~da/dt 
as a function of 7. The agreement with the decay 
rate computed using Chandrasekhar's formula ^TS\\ (open 
squares) is good. For 7 = 0.6, there is not any significant 
evolution of the orbit in the considered interval of time 
and, consequently, s « 0. 

A similar conclusion is implied by Figure |9] which shows 
the trajectory of a 2OOOM0 black hole in model D, a 
longer integration with A^ — 130000 and 7 = 0.6. Ini- 
tially, the black hole sinks rapidly to the center, reaching 
~ Tst in ~ 3Myr. As the inspiral progresses, the orbit be- 
comes more eccentric (e « 0.3 at 4Myr). At later times 
(> 4Myr), the orbit shows no sign of further decay, oscil- 
lating in radius between ~ 0.1 and ~ 0.2pc. The orbital 
eccentricity remains almost constant in this phase. 

These findings, obtained for a flattened density cusp 
around a SMBH, seem to confirm the theoretical pre- 
dictions made above: i) dynamical friction "vanishes" 
within ~ O.lSpc; ii) the orbital eccentricity of an in- 
falling body increases with time. 

However, in any A'^-body simulation, stars are contin- 
uously scattered by gravitational encounters with other 
stars, with the result that the initially empty phase space 
region responsible for the vanishing dynamical friction 
force will gradually be filled. In addition, due in part 
to the low central density of our GC models when 7 
is small, the radius at which the cumulative mass in 
stars becomes comparable to that of the inspiraling black 
hole can be of order rgt, even for relatively small M 
(see table 1). A^-body simulations have shown that, in 
these circumstances, the orbit deviates from the theo- 
retical prediction of the Chandrasekhar's formula as a 
consequence of perturbations induced by the infalling 
black hole on the inner cusp (iBaumgardt et al. 1 120061 : 
iLockmann fc Baumgardt"1l2008[ ). Finally, it is not clear 
whether the approximations made in deriving equation 
(fT8|). which was the basis for the red lines plotted in 
Figure m are reasonable, or how large might be the fric- 
tional force from fast moving stars that populate the low 
density core. In fact, as we now demonstrate, these ad- 
ditional effects have a substantial infiuence on the long- 
term evolution of the black hole orbit. 
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Fig. 8. — Left panel; evolution of the semi-major axis for a 5000Mq black hole in the short A'^-body integrations, for different values of 
the central density slope (from top to bottom, 7 = 0.6, 1, 1.5, 1.8). Th e th icker line is from the high-A^ integration, with A' = 500,000 
and 7 = 0.6. Dashed lines are predictions from Chandrasekhar's formula l|18|l using InA = 6.6. For 7 = 0.6 there is no significant evolution 
of the orbit in the considered interval of time. Right panel: orbital inspiral rates s = —da/dt computed for the simulations displayed on 
the left panel as a function of 7 (filled circles). Open squares give the predictions from Chandrasekhar's formula. The star symbol is the 
decay rate computed from the high resolution run {N = 500, 000 and 7 = 0.6). 
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Fig. 9. — Trajectory of a 2000Mq black hole into a core vfith 
7 = 0.6 (model D). The top-red line is the theoretical prediction 
obtained from Chandrasekhar's formula l|18|l using InA = 6.6. The 
bottom red curve shows the predicted inspiral in a 7 = 1.8 cusp. 



Figure [TU] shows the trajectory of the black hole 
for some of the iV-body integrations from Table 1 and 
compares them to the evolution predicted by Chan- 
drasekhar's formula ()18|) (upper green curves). (In the 
upper panels, the comparison is displayed only for the 
higher resolution runs, i.e., models Al and Bl.) Al- 
though the agreement with the theoretical prediction 
appears fairly good, at least for M = SOOOMq, when 
7 = 0.6, the A^-body integrations reveal a faster decay 
than predicted. Either some the frictional force must 
come from stars with velocities > w, or the back- 
ground stellar distribution is changing during the inspi- 
ral (or both). These two possibilities are investigated in 
what follows. 

Dynamical friction from fast-moving stars. Equation 
(fT8|) was derived under standard approximations that ig- 



nore the contribution from non-dominant terms and the 
velocity dependence of InA. Although these approxima- 
tions are reasonable when there is a large fraction of stars 
with low velocities (i.e., v^, < v), it is unclear whether 
they can be applied to a region populated mostly by stars 
moving faster than the black hole. 

Without these assumptions, the instantaneou s dy- 
nami cal friction acceleration becomes (jChandrasekhar I 
[Toil : 



ff,^~4TTG^Mp{r)— / dvATTf{v,)vl 
Jo 



\v — v^ 



dV 1 



In 1 



max * 



(29) 



where /(w*) is the velocity distribution of field stars, 
and Prnax IS the effective, maximum value of the im- 
pact parameter. In this more accurate treatment, some 
of the dynamical friction force is due to stars moving 
more rapidly than t h e massive parti cle (jChandrasekharl 
[T94I IWhite I [T949I: IMerritt I [200l . If the condition 
PmaxV^ /GM » 1 is satisfied, the frictional force can be 
approximated as (Chandrasekhar 1943, equation 30): 



ffr' 



' •/ fr "T •/ fr 



-47rG^Mp(r)- 



dvA'^f{v*)vl In 



GM 



(30) 



dvA'^f{Vi,)vl 
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Inside rst , dynamical friction is produced mostly by stars 
with Vi, > V and the first term in the integral becomes 
negligible. This shows the weak dependence of the fric- 
tional deceleration inside the core on Pmax- 




cio 5x10^ 10^ cio 2x10* 4x10^ e-KlO^ ej<10^ 10'' 



Tinne(yr) 




Tlme(yr) 

Fig. 10. — Orbital evolution of the second black hole in models Al, A2, Bl, B2, E and F. Solid green lines show predictio ns a ssuming 
a fixed background of stars. Upper green curves are obtained by us ing the standard Chandrasekhar's formula (i.e., equation l|18|l ). while 
lower green curves give the orbital decay computed using equation II29I I with pmax = 0.5 pc. Red lines were obtained with equation II29I I 
but allowing /(t)*) and p{r) to change according to the evolution of the A^-body system. 



Adopting equation (P^ . with Pmax — 0.5pc, for the 
frictional force that appears in the equations of motion 
(|20l). we obtained the lower green curves in Figure fTOl 
which show much better agreement with the A^-body re- 
suhs. Evidently, the standard expression for dynamical 
friction, equation (jlSp , is inadequate to describe the or- 
bital evolution of a massive body at the GC in the case 
that the density profile of the nuclear star cluster is shal- 
low. This is apparently a consequence of neglecting the 
non-dominant terms, and not, for instance, of the as- 
sumed independence of the Coulomb logarithm on the 
field star velocity distribution. For models Al and A2, 
Lagrangian radii showed essentially no evolution, indi- 
cating the absence of any significant change in the stellar 
distribution induced by the second black hole. We con- 
clude that (at least) some of the drag within Tgt is due to 
field stars with Vi, > v. The red lines in Figure [TOl were 
derived from equation (j29p but using a time dependent 
distribution function f{vi,,t) extracted (at time t) from 
the A^-body models (see below). For models Al and A2 
the red curves agree exceptionally well with the A^-body 
results and they essentially match the results of the semi- 
analytical integration that takes into account the friction 



from fast moving stars. We conclude that for these runs 
it would be appropriate to ignore the influence of the 
second black hole on the stellar distribution. 

In the left panel of Figure [TT] we plot the fraction of 
the dynamical friction force that is predicted, by equa- 
tion (l29t . to come from stars with v^, > v, for different 
values of the inner cusp slope and at different radii. In 
the right panel of the figure, the total frictional decel- 
eration in our models is given in units of the frictional 
force computed under the assumption of a Maxwellian 
distribution of velocities: 



/fr 



-47rG2Mp(r)lnA 



erf(A) - —e 



(31) 



with X = v/\/2a. Clearly, this equation, often used in 
the past to describe the orbital evolution of a massive ob- 
ject into the GC, overestimates the frictional drag within 
r < 0.2pc for 7 < 1. 

Influence of the second black hole on the field-star dis- 
tribution. For larger masses of the infalling body, i.e. 
Al > IOOOOM0, the perturbations which it induces in 
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Fig. 11. — Left panel: fraction of the dynamical friction force that is predicted to come from stars with t;*^ > t; as a function of 7, at 
different galactocentric radii: r = 0.1,0.2,0.3 and 0.6 pc. Equation II29I I was used to compute these curves. When 7 = 0.6, dynamical 
friction at small radii comes only from stars with 11* > v. As either 7 or r increase, the contribution from fast moving stars decreases. 
Right panel: total dynamical friction force in units of the frictional deceleration computed assuming a Maxwellian distribution of velocities. 
The frictional force produced by stars with Vi, > v, in the flattened cusp (i.e., 7 = 0.6 and r < 0.2pc) is much smaller than that obtained 
under the simple assumption of thermal distribution of velocities. In both panels we adopted Pniax = 0.5pc and M = IOOOA/q. In the right 
panel, we used InA = 6.6 to solve equation |(3TJ. 
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Fig. 12. — Left panel: Fraction of stars with velocities less than the local circular velocity F{< tJcirc^") as a function of radius, at the 
same time (3 X lO^yr) for models A2 (M = 5OOOM0), E (lOOOOAf©) and F (5OOOOM0) . The dashed curve corresponds to the initial 
configuration. The larger the mass of the black hole the faster the changes of the model in velocity space. Right panel: F{< fcirc,^) as a 
function of radius for model Al at different times. Due to two-body relaxation, stars are scattered toward low velocities and the hole in 
phase space that characterized the initial configuration is gradually filled up. 



the background system introduce a complex time depen- 
dence of the phase-space distribution. During the orbital 
inspiral, the black hole scatters stars into the inner cusp; 
consequently, once it reaches ~ Vst, it will "see" stars 
with Vi, < V that contribute to the frictional acceleration 
from that point on. 

In order to test Chandrasekhar's formulae under these 
circumstances, the black hole equations of motion were 
integrated in a time-varying potential whose properties 
were varied over time in a way designed to mimic the 
evolving iV-body models. In more detail, the density of 
the A^-body model was computed at fixed intervals of 



time by binning particles in concentric logarithmically- 
spaced shells. At the same time the velocity distribu- 
tion of field stars was obtained directly from the A^-body 
model. Finally, the black hole equations of motion were 
numerically integrated as described in section (|3.ip using 
expression (P^ . In this way, we were able to approxi- 
mately account for the back-reaction of the second black 
hole on the stellar distribution. It is worth noting that, 
even with this more sophisticated approach, two relevant 
assumptions are retained: i) any induced deviation of the 
models from isotropy is neglected; ii) the black hole is as- 
sumed to move always on a circular orbit, while the A''- 
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Fig. 13. — Lagrangian radii evolution of models F (upper panel) 
and E (lower panel). Green curves show the position of the massive 
body. 

body simulations clearly show an increase of the orbital 
eccentricity with time. The red curves of Figure [TOl 
obtained through this numerical procedure, show that 
even when the galactic nucleus is rapidly deviating from 
its initial configuration, Chandrasekhar's theory can still 
accurately reproduce the TV-body results if the changes 
in the stellar distribution are taken into account and the 
fast moving stars are included when computing the fric- 
tional force. 

In Figure [T^lwe show the evolution induced by the sec- 
ond black hole in the velocity distribution of the model, 
by plotting the function F{< Wcirc,?') at the same time 
(3 x lO^yr) for different masses (left panel). In addition, 
we show how F(< Vcirc,r), for M = SOOOM©, evolves 
as a function of time (right panel). In this latter case 
two-body relaxation causes the diffusion of stars at low 
velocities and the stalling radius is shifted from the ini- 
tial ~ O.lpc to ss O.OSpc by the end of the simulation. 
We note that - in a real galaxy with much larger N - 
this effect would be essentially absent. 

Figure [13] illustrates the changes in the configuration- 



space density for models E and F via the time evolution 
of their Lagrange radii. The time evolution of models 
E and F is remarkable: in model F, the perturbations 
on the stellar distribution are initially so large that the 
core fills up during the first ~ 2 x lO^yr. At this point, 
the black hole, at a galactocentric distance of ~ O.OSpc, 
starts to carve out the inner region, destroying the cusp 
that it created before. The final model has a core of 
size ~ 0.2pc and the internal slope is 7 < 0.5. However 
its density is, everywhere within Ipc, smaller than that 
of the initial model as a consequence of displacement of 
stars from the cusp. A qualitatively similar evolution was 
found in model E. Figure [HI shows the induced evolution 
of the density profile for runs E and F as well as the time 
variation of the anisotropy parameter, defined as 

/3 = 1 - a^/a^ , (32) 

with at and cr^ tangential and radial velocity dispersions 
respectively. 

In summary, a straightforward interpretation of our A^- 
body results is that equation ([T^ reproduces remarkably 
well the real decay rate of a massive object into the GC 
only until it reaches the stalling radius. In the subsequent 
evolution, the orbital decay slows down as a consequence 
of the lack of slow moving stars in the inner galactic 
nucleus (see Figures [8] and [9]) , but it never drops to zero, 
due apparently to the frictional force generated by stars 
moving faster than the inspiraling black hole ( Figures [TU] 
and[II]). 

A massive body of mass M 1000 A/©, starting from 
distances of order rbh, will reach a galactocentric ra- 
dius ^ O.Olpc in ^ lO^yr. For larger masses (i.e., 
M > IOOOOA/0), during the inspiral, the black hole en- 
hances the diffusion of stars into the phase-space region 
that was initially nearly empty (Figures [12] and [T3l) . Dur- 
ing the stalling phase a low density core is rapidly re- 
generated by the second black hole as it displaced stars 
from the cusp. Notice that, in our models the stalling 
distance is about ten times larger than that found in pre- 
vious works that assumed a collisionally-relaxed, steeply- 
risin g density profile around the central black hole 
(e.g.-lBaumgardt et al. 1120061 iLockmann fc Baumgardt I 

mm ). 

We note in passing that the background stars have 
orbital periods similar to that of the massive body. It 
is conceivable that correlations may be induced by the 
massive body in the orbital elements of the stars that 
will change the evolution significantly away from that 
produced by an uncorrelated background. On the other 
hand, two-body relaxation in the A^-body models will 
tend to de-correlate the background response, leading, 
perhaps, to a better correspondence with the predictions 
of Chandrasekhar's theory. 

4.3.2. Eccentric Orbits 

In this section, we investigate the rate of change of the 
orbital eccentricity as a consequence of dynamical fric- 
tion. We devised two simulations that differ only in the 
number of particles: 200,000 and 100,000. We refer to 
these simulations as runs Gl and G2 respectively (see 
Table 1); both have 7 = 0.6. The black hole was ini- 
tially placed at a radius of r-m = Ipc on an eccentric 
orbit with Cm = 0.54. As discussed earlier (Section [3?2]) . 




log r(pc) log r(pc) 

Fig. 14. — Left panels: density profile evolution in run F (upper panel) and E (lower panel). The black curve corresponds to the initial 
model; the red line is obtained at time lO'^yr for run E and at 2 X lO'^yr for run F, while the blue lines are the density profile of the final 
models, after the secondary black hole has stalled carving out a deficiency of stars in the inner regions. Filled circles indicate the position 
of the inspiraling . Right panels: Evolution of the anisotropy parameter in the models. Line thickness increases with time. As the black 
hole spirals in, it induces tangential anisotropy in the background system. 



when the orbital periapsis hes within the core, the orbit 
is expected to become more eccentric as a consequence 
of the declining frictional force in this region. 

Figure 1151 shows the evohition of the eccentricity and 
semi-major axis of the orbit as a function of time, demon- 
strating that, at least quahtatively, Chandrasekhar's the- 
ory reproduces the evolution. Although the eccentricity 
undergoes significant fluctuations, it evidently drifts to- 
ward larger values with time. This behavior is quite ro- 
bust showing a negligible iV-dependence. 

It is generally assumed that dynamical friction, in 
power-law density models with an isotropic velocity dis- 
tribution, would circularize the orb i t of a n infalling body 
(see for instance lBaumgardt et al. I ()2006D ). Our A^-body 
simulations demonstrate that in models characterized by 
a flat density profile and a central SMBH, the eccentricity 
can instead be an increasing function of time. 

5. GRAVITATIONAL WAKE 

An alternate way to look at dynamical friction is 
in terms of the acceleration produced by the over- 



density of stars that accumula te behind the massive 
body - the "gravitational wake" ( Danbv fc Camm lllQSTt 
iMarochnik] Il968l: iMulder 1 119831 1 The expression for 
the response wake in a homogeneous mediu m is given 
for arb itrary spherical density distribution in iWeinbergI 
()1986[ ). The existence of a wake has rarely been 
confirmed in A^-bo dy simulations; an isolated exam- 
ple is provided by iWeinberg fc Katz I ()2002D (see also 
[Weinberg fclcitil (I2007D 1 who show the wake induced in 
a dark-matter halo by a stellar bar. Other examples in- 
clude Weinberg (1989} , Hcrnquist & Weinberg (1989), 
iVesperini Weinberg I (j2000l ). 

We searched for the wake in our A^-body simulations by 
computing the relative overdensity at each radius along 
the orbit of the second black hole. The A^-body models 
were first rotated in such a way that the second black 
hole was situated at y = z = with = and Vy > 0. 
The density at any position was then estimated using a 
Gaussian kernel with radially- varying smoothing length. 
Figure [12] shows the results in runs Al, E and F as a 
function of the azimuthal angle 9 at different radii and 
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Fig. 15. — Evolution of eccentricity and semi-major axis for models Gl (upper panel) and G2 (lower panel) that differ only in the number 
of fi eld particles: A'=200,000 and 100,000 for mod els Gl and G2 respectively. Dashed lines are the theoretical predictions from equation 
ITSt . Dotted lines were obtained with equation 11291 1 (i.e., including the frictional drag from stars with > t?), where we used pmax = O.Spc. 
As the black hole spirals in, its orbital eccentricity increases. This conclusion is quite robust, showing essentially no dependence on the 
number of background particles. 



for different values of M. In the figure, the black hole lies 
at 9 — with 6 > and the average density is defined 
as (l/27r) J^^dO p{9) : outside the core (r > 0.3 pc), 

the peak in the overdensity lies at —20 < 9 < 0°, in- 
dependent of Af, and the amplitude of the overdensity 
increases with black hole mass, as expected. The wake 
is therefore always just behind the massive body in this 
phase. When r < 0.3 pc, for M = 5000 - IOOOOMq, 
the density enhancement is reduced but its position re- 
mains essentially unchanged. The reduced amplitude of 
the wake inside the core explains why the frictional force 
is greatly suppressed in these regions. For larger masses, 
the angular dependence of the overdensity in this phase 
is more complex, revealing, in some cases, two distinct 
peaks. During this phase, the mass distribution is af- 
fected by gravitational scattering from the massive body. 
Finally, when the black hole is well inside the core, the 
density maximum is seen to lie at large angular separa- 
tions {9 < —100°) from the black hole. Indeed, a den- 
sity "hole" , with amplitude approximately proportional 
to M, is apparently induced by the black hole at roughly 
its position during the stalling phase. 

Figure [T7l shows two-dimensional contour maps of the 
overdensity for run E (M = 5 x 10**). The radial ex- 
tension of the wake (with respect to the galaxy center) 
does not change greatly over time, but one can clearly 
see how the location of the density maximum shifts, and 
a density gap is apparently created near the black hole 
position during the stalling phase. 



To more clearly illustrate how the location of the grav- 
itational wake with respect to the second black hole 
evolves, we plot in Fig [TH] the angular position of the 
maximum as a function of the black hole galactocentric 
radius. Outside the core (i.e., r > 0.3 pc) the wake is lo- 
cated at small (negative) angles, causing the initial rapid 
inspiral. Once the black hole starts to modify the back- 
ground of stars the wake becomes more difhcult to track. 
This causes the large oscillations seen in the relative po- 
sition of the wake and in turn explains why such oscil- 
lations occur earlier for larger masses of the inspiraling 
object. 

6. DISCUSSION 

In this paper, we presented TV-body simulations of 
the inspiral of a massive body into the Galactic center 
(GC). Our models of the Milky Way nuclear star clus- 
ter were motivated by recent observations that suggest 
a relatively low density of stars inside the SMBH influ- 
ence radius. Such models are characterized by a zero 
or near-zero phase-space density at low energies. Un- 
der the standard approximation, in which the frictional 
force from fast-moving stars is ignored, a second black 
hole that sinks toward the center under the influence of 
dynamical friction would stall at a distance of roughly 
1/2 the core radius, or ^ 0.25pc, from the SMBH. If the 
smaller black hole moves initially on a non-circular orbit, 
its orbital eccentricity is predicted to increase with time 
due to the lower dynamical friction force near periapsis. 

Using TV-body simulations, we found that the fric- 
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Fig. 17. — The density response (i.e., gravitational wake) induced by the massive body in run E is shown by plotting density contour 
maps of background stars in the upper panels, and the corresponding relative overdensity along the black hole orbit in the bottom panels. 
The isodensity contours were obtained by subtracting at any radius the mean density and selecting only particles that were close to the 
orbital plane. Negative contours (underdensities) are shown by dashed curves. Circular regions show the path over which the density was 
computed to obtain the plots in the bottom panels. 



tional force never falls precisely to zero. As noted also 
by Chandrasekhar, stars moving faster than the test 
body contribute to the drag. When this contribution 
is included in the expression for the dynamical fric- 
tion, Chandrasekhar's formula reproduces quite well the 
decay observed in TV-body simulations of the inspiral 
of a ~ IOOOM0 black hole. The eccentricity increase 
predicted by Chandrasekhar's theory is also confirmed. 
When the inspiralling object is more massive, a second 
mechanism contributes to the frictional force: the second 
black hole induces evolution of the background system, 
which tends to refill the initially empty regions of phase 
space. 

In what follows, we discuss the implications of these 
results for a number of astrophysical problems related to 
the dynamics of massive bodies near the centers of galax- 



ies. But first, we comment on how our iV-body results 
can be approximately scaled to systems with different 
masses and densities. 

The rate of inspiral of a massive body of mass M is 
independent of the mass of field stars if M ^ m. Chan- 
drasekhar's formula alsopredicts a linear dependence of 
the frictional force on AiQ and our simulations (as well as 
many others) confirm that prediction. If the density re- 
sponse of the background is ignored, the iV-body results 



^ In its more general form JTSj, the dynamical friction formula 
predicts an additional, approximately logarithmic dependence of 
force on M . 
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Fig. 18. — Position of the relative density maximum as a function 
of the black hole galactocentric radius in runs Al (fille d circles), E 
(open circles) and F (stars symbols). As in Figure [161 the A'^-body 
models were rotated such that the second black hole is located at 
8 = with e > 0. 
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where r is radius containing a mass in stars Af^,(< r) « 
3.4 X 10~^M,, M is the mass of the test body and M 
its mass adopted in the iV-body simulations of Table 1 . 
When in the simulations the background of stars evolves, 
the dependence on the mass of the infalling body be- 
comes more complex; in this case, the appropriate scaling 
is obtained by setting M = M x [M,/A x 10*^ Mq] , i.e., 
setting the ratio between the mass of the massive body 
and the central black hole the same as in the A'^-body 
simulation.^ Particular caution should also be taken when 
adopting M > M since for large values of M the mas- 
sive body would perturb the background system more 
importantly than it does in A'^-body runs. 

The condition that the background not evolve is satis- 
fied in our simulations when M < IO'^A/q and at early 
times in run E. We apply this approximate scaling to run 
Al, for which M = SOOOM© and the total integration 
time is ~ 1.5 x 10^ yr. Assuming no change in the stellar 
density, replacing the massive body by a ~ 10 Mq black 
hole increases the effective integration time by a factor 
^ 500, to ~ 8 x 10^ yr (at which time the galactocentric 
radius is ~ 0.06 pc). This result illustrates how - in the 
absence of a steep central density cusp ~ the time for 
stellar-mass BHs to reach the center of the Galaxy from 
a starting radius of ~ 1 pc can easily exceed ~ lOGyr (a 
point we return to in Section [6.ip . 

Alternatively, we can identify our models with the cen- 
ter of a galaxy like M87, a luminous elliptical galaxy with 
a flat central density profile. We adopt M, = 3 x 10^ Mq 
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Fig. 19. — Evolution of the density profile of a population of 
10 Mq BHs (dotted curves) assuming that they constitute 1% of 
the total mass density initially. Results are displayed for three 
choices of the core parameter ro =(0.3, 1, 2) pc. Lower (upper) 
solid lines show the initial density profile of stellar BHs (stars). 
In the upper-left panel the BHs lie on circular orbits while in the 
other cases we assume an isotropic initial distribution of velocities. 
Density profiles are shown at time intervals of At = 2x10^ Gyr in 
the lower panels, while At = lO^Gyr in the upper-right panel. 
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ity dispersion cr„ = 278kms^^ and the relation ct,^ 
47rG'/9o(ro/3 )^ with rp = 600 pc to obtain the mass den- 
sity profile ([Young et al. II1978I: iLauer et al. II1992D : 



p{r) = 35 Mopc 



-3 ( 



V600 pc 



(35) 



Taking 7 ~ 0.6, this gives a length normalization factor 
f « 20 pc. 

Runs E and F: in these runs the background sys- 
tem evolves due to the perturbations induced by the 
massive body (see Figure 14). Setting M = M x 
[M,/4 X lO^M©] , run E corresponds to the inspiral of 
a 7 X 10^ Mq black hole starting from a distance of 
200 pc, and a total integration time '--^ 2 x lO^yr. In the 
case of run F, the inspiraling black hole would have a 
mass ~ 4 X 10'' M©; it penetrates the inner ^ 10 pc in 
~ 3 X 10* yr after which it effectively stalls. 
Run A: the condition that the background not evolve is 
satisfied in runs Al, A2 and also at early times in run 
E. Setting M = 10^(105)10-* Mq in runs Al and A2, 
the final integration time and orbital radius are ^ 3 x 
10®(10^'^)10^^ yr and 12 pc respectively. This shows how, 
in the central core of a M87-like galaxy, the inspiral time 
for black holes of masses 10^ Mq could easily exceed a 
Hubble time (a point we further discuss in Section !!). 2. 2p . 

6.1. Segregation of massive remnants at the Galactic 

center 

About 1% of the total mass of the old population 
at the GC should be in the form of stellar-mass (to ~ 
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Fig. 20. — Left panels: final cumulative eccentricity distribu- 
tion of stellar BHs for the integrations displayed in Figure 1191 (red 
curves), that would be measured inside the core within some radius. 
Solid curves give the initial distributions. Right panels: cumula- 
tive eccentricity distributions of the initial models (solid curves) 
evaluated within different radii. At small galactocentric radii, the 
distribution is dominated by high eccentricity orbits, in spite of the 
fact that the velocity distribution is isotropic. Dashed curves show 
for comparison a "thermal" eccentricity distribution, A' ~ e^. 



10 - 20Mq) BHs () Alexander I [2005[ ) . Since stellar BHs 
are significantly more massive than the mean stellar mass 
(~ IMq) expected for an evolved population, they would 
spiral in to th e center and segregate around the SMBH 
(jMorris Ill993fl . The time evolution of the remnant pop- 
ulation depends sensitively on its initial distribution and 
also on the properties of the background distribution of 
lighter stars. In the case of a flat core in the stars, and 
a similar initial distribution in the BHs, the time for 
the latter to reach a steady state density profile can ex- 
ceed a Hubble time, since the dy namical friction force 
essentially ceases inside the core (IMerrittllMol) . On 
the other hand, if the stars follow a steep central den- 
sity cusp, the mass density of BHs after ^ 10 Gyr can 
reach or exceed that of the other populations within 
^ 10~^ pc, leading to a qua si-steady-state density profile 
n ~ at small radii fe.g.lHopman &: Alexander"1l2006l : 
[Alexander k. Hopman 1120091 1 

Understanding the distribution of BHs at the centers 
of galaxies like the Milky Way is crucial for making pre- 
dictions about the expected e vent rate for lo w- frequency 
gravitational wave detectors (jHughes |[2003l ). Since the 
stellar BHs at the GC are not directly detected, time- 
dependent inspiral calculations like the ones presented 
here provide the best hope of understanding their distri- 



bution. However, if the background stellar distribution 
is a flat core, our results show that a straightforward ap- 
plication of Chandrasekhar's formula can give misleading 
results. 

Accordingly, we computed the evolution of a popula- 
tion of stellar BHs as they spiralled in to the center of 
a galaxy with a fiat stellar core, including the frictional 
force from the fast-moving stars. We began by generat- 
ing random samples of positions and velocities from the 
isotropic distribution function corresponding to the den- 
sity model of equation (jl4p assuming 7 = 0.6; cores of 
various sizes, ro — (0.3,1,2) pc; and selecting only par- 
ticles within 5 pc of the SMBH. In each of these models, 
a total of 800 orbits (representing the stellar BHs) were 
then integrated by solving the system of equations (|20)) . 
with dynamical friction force given by 



Si^ = -^'kG^M p{r)—[F{< v,r) In A 



(36) 



dvA'^f{Vi,)vl 



In 



V* - V 



with InA = 15, AI — 10 Mq. At each time, the density 
profile and eccentricity distribution of the inspiralling ob- 
jects were computed by sampling each orbit over time 
intervals of 0.3 Gyr. We also considered one model with 
core parameter tq = 0.3 pc in which all BHs were initially 
on circular orbits. 

All of the calculations presented in this section assume 
that the mass density due to the BHs remains small com- 
pared with the mass density in stars, and that the stellar 
distribution is unchanging. Because the two-body relax- 
ation time for 1 Mq stars is so long in these models, 
and ^10 times longer than the black hole inspiral time, 
ignoring the evolution of the stellar distribution due to 
star-star encounters is reasonable. This basic assumption 
is also supported by recently publis hed A^-body simula- 
tions ([Gualandris fc Merritt 1 12011[ ) that show how, in 
models with a pre-existing stellar core, the distribution 
of BHs evolves against an essentially fixed background of 
stars. However, once the density in BHs begins to ap- 
proach that in the stars, our calculations are no longer 
valid. 

In Figure[T2]we plot the density profile of BHs at differ- 
ent times, assuming that their fraction is initially 10~^ 
of the total mass density. The upper panels give the 
results for the model with rg = 0.3 pc. In these inte- 
grations the core is very small and after only ^ 1 Gyr 
the density of black hole rises very steeply going into the 
stellar core. After ~ 4 Gyr the BHs accumulate at radii 
near the core, matching the density in stars at ~ 0.01 pc. 
In the circular-orbit model, the density profile at 1 Gyr 
shows a maximum at ^ 0.2 pc, that grows and migrate 
inward due to the friction produced by fast moving stars 
inside these radii. The evolution for the isotropic run is 
comparably rapid, and after ^ 3 Gyr the density of BHs 
r eaches that in s tars at ~ 0.01 pc. 

IMerritt I (j2010 D showed that a core of the size currently 
observed is a natural consequence of two-body relaxation 
acting over 10 Gyr, starting from a core of radius ~ 1 pc. 
It is therefore of interest to study the evolution of the 
black hole distribution in models with parsec-scale cores. 
This is shown in the lower panels of Figure [191 In these 
cases the evolution is slower as a consequence of the in- 
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creased size of the stellar core, and even after 10 Gyr the 
density of BHs can remain substantially lower than that 
in stars at all radii. We conclude that it would be un- 
justified to assume that the massive remnants have yet 
reached a steady-state density profile at the GC. One 
consequence is that rates of capture of stellar BHs by 
the supermassive black hole at the Galactic center (EM- 
RIs) may be much lower than in standard models that 
postulate a collisionally-relaxed nucleus (e.g., Hopman & 
Alexander 2006). 

The left panels of Figure [20] plot the cumulative distri- 
bution of eccentricities of BHs inside various radii. Since 
the final eccentricity of each orbit is larger than its ini- 
tial value (see Section 2.2), one might naively expect the 
eccentricity distributions to evolve toward a form that 
is increasingly strongly peaked near e « 1. This would 
be the case if one plotted A^(e) for a fixed subset of ob- 
jects. However, when restricting the sample to a given 
radial range, t he result is very different. The reason (e.g. 
iMerritt 1120101 Appendix) is illustrated in the right-hand 
panels of Figure [20l given a flat density proflle, even 
an isotropic distribution of objects around a SMBH will 
have an eccentricity distribution that is strongly peaked 
near e = 1, since the only objects that can approach 
closely to the SMBH are on highly eccentric orbits. As 
the distribution of BHs evolves away from this initial 
configuration, the regions of low-energy phase space that 
were initially empty are gradually refilled, and the ec- 
centricity distribution begins to approach more closely 
to a "thermal" form, iV(< e) oc e^. In addition, (i) 
the eccentricity of individual orbits inside the core grows 
only very slowly since they are in a region where the 
dynamical friction force is small (see Figure 3); (ii) the 
eccentricity of BHs initially beyond the core decreases 
initially since they lie in a 7 « 1.8 cusp; their eccentrici- 
ties subsequently increase as the orbital periapsis enters 
the core, but in most cases this second phase is too short 
(see Figure 4) to produce final eccentricities significantly 
different from the initial values. We finally computed the 
anisotropy parameter (|32|) at the final integration time, 
defined as the time when the mass density in BHs reaches 
that in stars at small radii, and found that the departures 
from isotropy remained small at all radii. 



6.2. Dynamical Evolution of Eccentric Black Hole 
Binaries 

Gravitational radiation emitted by binary black holes 
with masses 10"^ — 1O^M0 is the principal target of 
planned, space-based, gravitational wave observatories. 
In the present literature the strain amplitude of the gravi- 
tational wave (GW) radiation is typically obtained under 
the assumption of complete circularization of the binary 
at the moment that the signal enters into the observ- 
able band. This simplification is motivated by the pre- 
dicted strong eccentricity decay when binary dynamics 
are dominated by relativistic effects. The expressions of 
the time average change of eccentricity e and semi-major 
axis a in the relativistic regime of a binary with com- 
ponents of masses toi and m2 were derived by Peters & 



Matthews (1963): 

/ da\ 64 G"^ TOi TO2(toi -|- TO2) 



dt/ 5 c^a3 /(e), (37) 

de\ _ 304 G^mim2{mi+m2) / , 121 3\ /oon 
dt/ 15 c5a4(i_e2)5/2 [^+304^)^ ' 

where c is the speed of light and 

/(e) = (l-e^)-^'^(l + ge^ + ge^) . (39) 

The strong dependence of the enhancement factor /(e) 
on e, shows the fundamental role of the binary eccen- 
tricity in determining the rate at which the system loses 
energy due to GW emission. 

A way to follow the orbital inspiral of a massive body at 
the GC, due both to dynamical friction and gravitational 
wave radiation, is to couple Chandrasekhar's formula for 
the frictional drag with the 2. 5 post-Newton ian equations 
representing GW energy loss (jMerritt Il2012t ) . In the limit 
M/M, << 1. the total deceleration can be approximated 
by 



f = -4:nG^Mp(r)^(F(< v,r) In A 



(40) 
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(41) 
(42) 



with Vj. the radial component of the velocity vector. Evi- 
dently, both the frictional force and the 2.5PN correction 
are dissipative terms, but, while the latter term always 
drives to lower eccentricities, the effect of dynamical fric- 
tion on the orbital eccentricity has a strong dependence 
on the phase-space distribution associated with the stel- 
lar background (see ii l3.2p . Notice that, in equation (|40| . 
if we neglect the dependence of the Coulomb logarithm 
on the mass of the test body, both the frictional term and 
the post-Newtonian terms depend linearly on M imply- 
ing that the time evolution of the orbital elements can 
be trivially rescaled to any M as long as the condition 
M/M, << 1 holds (see also equations 36 and 37). 

6.2.1. Dynamical Friction in the Context of the EMRI 
Problem 

Extreme-mass-ratio inspirals (EMRIs) are a poten- 
tial s ou rce of low-frequenc y gravi ta tional waves pughes 



2003t IBarack fc Cutler I I2004t lAmaro-Seoane et al. 



2007f ). In steady-state models of the Galactic center, the 
distributed mass with in 10^^ pc of the SMBH is d omi- 
nated by stellar BHs (jHopman fc Alexander! I200I . At 
these radii, dynamical friction is therefore typically ig- 
nored and it is assumed that captures for EMRIs are 
driven by gravitational scattering from other stellar BHs 
(e.g. iMerritt et al. , ,201 h) . On the other hand, if the 
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10"^ 10"^ 0.01 0.1 1 2x10^°4x10^°6x10^°8x10^° 

1 -e t (yr) 

Fig. 21. — Left panel: evolutionary tracks of a massive object in the Galactic center starting from various eccentricities = (0.3, 0.5, 0.7), 
from an initial semi-major axis = 0.2 pc and adopting two different inner slopes of the mass-density profile 7 = (1,0.6). Dot-dashed 
lines are the Schwarzschild barrier, equation II45I I. below which resonant relaxation is suppressed by relativistic precession and gravitational 
scattering is dominated by classical non-resonant relaxation. Vertical marks give the radii within which the two-body relaxation time 
scale for changes in angular momentum (4r,off ) becomes shorter than the time-scale of evolution for angular momentum in our integrations 
(*evol)i assuming 10 Mq for the mass of the inspiraling black hole. Inside these radii, for M < 10 Mq, our integrations are no longer valid 
since two-body scattering, rather than dynamical friction, would dominate the orbital evolution. For 6;^ < 0.5 and 7 = 1 (two right-most 
curves), at any radius, gff was always longer than tcvol B.nd no vertical marks are displayed. In the two left-most curves, the condition 
that tr eft > ievol radius wold instead require a slightly larger mass for the BH: M > 15 Mq. This shows that gravitational scattering 

from stars can be neglected and our integrations are valid for relatively small masses of the test particle. Within the Schwarzschild barrier, 
dynamical friction is therefore the main mechanism inducing creation of EMRIs. We also stress that in these integrations, changes in 
the stellar distribution are not taken into account. For instance, the stellar potential would be strongly perturbed when the mass of the 
inspiraling black hole becomes comparable to the mass in stars contained inside its orbital radius. As a reference, dotted lines in the panel 
display the radius within which the mass in stars in the model is 10 or 1000 Mq. Right panel: time evolution of periapsis (dashed lines) 
and apoapsis (continue line) for a 10 Mq BH. The sinking time scale decreases with increasing the initial eccentricity, and, for the set of 
computed orbits, it is shorter than 10-*^" yr only for £;„ = 0.7 (left-most curve in the panel). 



background stellar distribution has a flat core, the den- 
sity of BHs can remain small compared with the mass 
densi ty of other populations (e.g. Gualandris & Merritt 
I20T1I §6.1). Under these circumstances, at any radius, 
massive remnants might see a background whose density 
comes mostly from lighter stars and dynamical friction 
becomes a competing mechanism in driving capture of 
EMRIs. 

Using equation (j40|) we computed the trajectory of the 
test mass under a variety of assumptions for the back- 
ground system. Results of these integrations are dis- 
played in Figure [H] We considered orbits of different 
initial eccentricities ein ~ (0.3,0.5,0.7), starting from a 
semi-major axis — 0.2 pc. For the stellar background 
we used the density model of equation (14) with two dif- 
ferent values of the internal slope index: 7 = 1 (black 
lines) and 0.6 (blue lines). For an eccentric orbit in a 
flattened cusp, dynamical friction at apoapsis dominates 
the evolution causing a rapid increase of the orbital ec- 
centricity. In the simplified picture in which the fric- 
tional force at periapsis is vanishing small, the apoapsis 
distance remains unchanged in time, while the periapsis 
becomes progressively smaller; at some point, the mini- 
mum distance from the SMBH is small enough that the 
2.5PN terms start to dominate the evolution. The drag 
at periapsis then circularizes the orbit, and causes the 



merger of the two black holes. 

Near a SMBH, as long as the relativistic precession 
time scale is much longer than the orbital period, the 
mechanism that dominates the scattering of stars onto 
high-eccentricity orbits is resonant relaxation. Because 
in the potential of a point-mass the orbits are fixed el- 
lipses, perturbatio ns on a test particle a re not random 
but correlated (Ra nch fc Tremaine Ill996f ). The residual 
torque |T| « ^/NGm/r, exerted by the N randomly ori- 
ented, orbit-averaged mass distributions of the surround- 
ing stars, induces coherent changes in angular momen- 
tum AL = Tt on times t < tcoh, where the coherence 
time tcoh is fixed by the mechanism that most rapidly 
causes the orbits to precess (e.g, mass precession, rela- 
tivistic precession). The angular momentum relaxation 
time associated with resonant relaxation is: 

trr = f XT" — ) ^co/i , (43) 
V ^^coh J 

where Lc = ^/GM,a is the angular momentum of 
a circular orbit and |ALco/t| \Ttcoh\ is the accu- 
mulated change over the coherence time. Assuming 
that the precession is determined by the mean field of 
stars, the angular momentum relaxation time becomes 
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Fig. 22. — Eccentricity at the moment the binary enters the 
sensitivity window of planned space- based interferometers, Cf, as 
a function of the initial orbital eccentricity ei„ for the integra- 
tions displayed i n F igure 1211 Star symbols are for 7 = 1 (black 
cur ves in Figure [JTJ, empty circle for 7 = 0.6 (blue curve in Fig- 
ure [2TJ. The dot-dashed lines give Sf ignoring dynamical friction. 
For a given initial eccentricity and secon dary black hole mass we 
fixed the merger time by using equation 1 148 II and varying the ini- 
tial orbital semi-major axis. If we take a test mass of 10 Mq 



(1000 Mq) this corresponds to merger times of 10 , 10 



.10' yr 



(10^^, lO^'^, ...10^ yr) from bottom to top line. As comparison, the 
orbital eccentricity and merger time for the integrations of Figure 
1211 at the moment GW energy loss stars to dominate the evolution, 
are (from left to right of that figure) e ~ (0.9994, 0.998, 0.994, 0.97) 
and tm ~ (5.9 X lO'^, 1 X 10^,1.9 x 10^, 5 x 10^) yr. Horizontal 
line represents approximately the lowest value of that would 
require non circular templates for data analysis (e ~ 10"'*, Porter 
& Sesana 2010). 



(jRauch fc Tremaine |[T996[) : 



trr « 2.9X10V 



M, 



4 X lO6Af0 



1/2 



O.lpc 



3/2 



(44) 

Dot-dashed lines in the left panel of Figure [211 give the 
Schwarzschild barrier. Above these lines, resonant relax- 
ation is the most rapid mechanism affecting angular mo- 
menta; below the curves, relativistic precession becomes 
efficient at suppressing resonant relaxation and the grav- 
itational perturbations are dominated by classical "two- 
body" relaxation. The value of the angu lar momentum 
that defines the Schwarzschild barrier is (|Merritt et al. I 

unni): 
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e )SB ' 
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where N is the number of stars within radius a and Csb is 
a constant of order of unity. Beyond the barrier, the time 
for encounters to change the orbital angular momentum 



by of order itself is tr,eff = 2(1 — e)tr, where for the non- 
resonant re laxation time scale we adopt t he approximate 
expression (|Hopman fc Alexander ||2006[ ): 
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For our integrations to be viable, the time-scale for dy- 
namical friction to change the orbital angular momen- 
tum, iovoi ~ (1 — e)|d(l — e)/di|-^, must be shorter 
than ir,off at all radii. Vertical marks in the left panel 
of Figure [21] give the orbital radius within which tr.cff 
becomes smaller than tovoi assuming M — 10 Mq. For 
ein — (0.3, 0.5) and 7=1 (two right-most curves), at any 
radius, tr,off was always larger than t^voi E^nd no vertical 
marks are displayed. Increasing the initial eccentricity 
to Cin = 0.7, tr^cS cquals tevoi at ~ 0.03 pc. At smaller 
radii, two-body relaxation would dominate the orbital 
evolution and, for a 10 Mo black hole, this integration 
is not longer valid. Taking ein = 0.5 and 7 = 0.6 (blue 
curve in the figure) this transition occurs at ~ 0.02 pc. 
Because increasing the mass of the test body reduces 
iovoi but leaves tr^eS unchanged, it is always possible to 
set M such that the condition ir,off > ^cvoi is satisfied 
everywhere within the Schwarzschild barrier. In these 
two latter cases this condition requires a slightly larger 
mass of the test body: M > 15 Mq (Woosley e t al~l 
l2002f ). Two-body scattering effects from field stars can 
therefore be ignored for relatively small masses of the 
sinking black hole. We conclude that, in a flat density 
distribution near a SMBH and at radii relevant for the 
EMRI problem (< 0.01 pc), dynamical friction might be 
an important process in driving the formation of EMRIs. 

Gravitational scattering can be dominated by other 
stellar black holes if their density becomes compara- 
ble of that in stars at small radii as a consequence of 
mass-segregation. In an unsegregated model the number 
of stellar black holes (of mass 10 Mq) is predicted to 
be 10"'^ times that in stars. From equations (|45p and 
(f46|) it follows that, in this case, the scattering from 
black holes can be ignored with respect to the pertur- 
bations from the stellar population. Gravitational scat- 
tering from black holes starts to compete with that from 
stars when their number at small radii (~ 1 mpc) is 
10^^ X N, similar to the found at later times in Fig- 
ure for rg ^ 2pc. In relaxed mass-segregated models, 
instead, the number of black holes would be approxi- 
mately N, and they will dominate the orbital evolution 
of the t est mass at any radius inside the Schwarzschild 
barrier ([Alexander fc Hopman |[2009f l. 

Finally, we note that dynamical friction can be very 
inefficient if the mass of the inspiraling object becomes 
comparable to the mass in stars within its orbital radius. 
In the 7 = 1 cusp for a ^ 1000(10) Mo, this occurs at 
- 0.02(0.002) pc or at - 0.03(0.005) pc when 7 = 0.6. 
This suggests that the results of Figure [2T] may not ap- 
ply for large masses of the test body and for small initial 
eccentricities (< 0.3). Accurate A^-body simulations, in- 
cluding high-order post-Newtonian terms, should be used 
to better understand at which extends the conclusions 
made here can be applied. We reserve this study to a 
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future paper. 

In order for an extra-galactic source to be observable 
by proposed space-based interfero meters, it must have 
an o rbital frequency > 10~'*Hz (jAmaro-Seoane et al. I 
[2001 . or 

/ M \ 
a < flf EE 4 X lO^Vpc \ . (47) 

We explored whether the computed orbits would retain 
some degree of eccentricity by the time the binary enters 
the instrumental sensitivity window, by evaluating the 
eccentricity, e^, at the time at which the condition (|47|) 
is satisfied and comparing this value with the minimum 
eccentricity that would require non-circular templates for 
data analysis: e ~ 10^* (Porter & Sesana 2010). We note 
that strong sources (with high eccentricity) might be de- 
tectable at l ower frequencies (i.e., larger semi-major axis) 
(jAmaro-Seoane et al. II2007I1 . The use of equation (j47| is 
therefore a conservative one. 

Figure [22] plots e/ as a function of the initial eccen- 
tricity for the orbits displayed in Figure [2TJ In addition, 
we computed a set of orbits with different initial eccen- 
tricities by removing from equation (j40p the dynamical 
friction term. Each dot-dashed curve in the figur e corre- 
spond s to a fixed value for the coalescence time (jPeters I 
[l96l : 




Taking M = 10 Mq (1000 M©), this corresponds to 
i„ = IQi^ yr (1013, 10^2, ...10*^ yr) from 

the bottom to the top line respectively. It is evident 
that even for relatively low initial eccentricities and large 
merger times the binary will have a value of e/ signif- 
icantly different from zero. This study suggests that 
secondary black holes typically reach the GW radiation 
regime on wide orbits that are still very eccentric, and 
even after the semi-major axis has decreased to values 
small enough for detection by space-based interferome- 
ters, eccentricities can be large enough that the efficient 
analysis of gravitational wave signals would require the 
use o f eccentric templates (see aslo iBarack fc Cutler I 
[200l . 

6.2.2. Orbital decay in the cores of giant elliptical galaxies 

Until the discqyery o f a stellar c o re in t he Milky Way 
(iBuchholz et al. I [2001 iDo et al. I [2001 iBartko et aU 
120101 ) ■ the densitv was generally assumed to follow a 
steep power law, p ~ r~^, inside the influence radius 
of Sgr A*. The same assumption is still commonly 
made when modelling the so-called "power-law" galaxies: 
galaxies of low to moderate luminosity that also exhibit 
steeply-ris ing densi t ies ne ar the center (jGebhardt et al. I 
[1996: Fa ber et al. I 119971 ). Whether other power-law 
galaxies will turn out to harbor parsec-scale cores like the 
one in the Milky Way remains to be seen. But it has long 
been known that cores are ubiquitous in stellar spheroids 
brighter than ^ 1Q^^L(7) , whose influ ence radii can be re- 
solved (jFerrarese et al. 1 ,1994; .Lauer et al. I I1995D . Core 



sizes are observed to be of order the SMBH influence 
radius or somewhat greater, consistent with models in 
which the cores are pr oduced by the scouring effect 
of bin ary SMBHs (jMerri tt 2006; iGualandris fc Merritt I 

[mil). 

In this section, we use equation (j40|) to investigate the 
orbital evolution of a massive black hole that spirals in 
to the center of a giant elliptical galaxy with a core. We 
base our models on M87. The relevant properties of M87 
are summarized at the start of this section. Here we 
note that the core of M87 extends substantially beyond 
the SMBH influence radius: Vc/rhh ~ 600pc/200pc w 3. 
By comparison, the Milky Way has Tc ~ 0.3rbh- This 
difference may reflect different formation processes for 
the two cores, or may be a result of the shorter relaxation 
time at the center of the Milky Way, which c ould cause 
the core to shrink over 10 Gyr (jMerritt 1 12010). 

Following the evolution of a binary SMBH at the 
center of a galaxy requires self-consistent simulations 
that can correctly treat the response of the back- 
grou nd stars to the pr esence of the second massive body 
(e.g. 'Ebisuzaki et al. "1991'; ' Quinlan fc Hernguist lll997l: 
Milosavljcvic fc Merritt 200l|). Here, we limit ourselves 
to the case where the inspiralling black hole is much less 
massive than the central one. For instance, capture of 
a Milky- Way-sized galaxy by M87 would bring a sec- 
ond SMBH into the center forming a binary of mass 
ratio ~ 10^"^. This problem can be seen as a scaled- 
down version of the capture of an intermediate-mass 
black hole by Sgr A*. Simula tions of the latter scenario 
(e.g. iBaumgardt et al.1l2006[ ) have generally assumed a 
steeply-rising stellar density around the SMBH; inspiral 
of the intermediate mass black hole is found to stall when 
the semi-major axis of the binary drops to ~ 10"'^ pc, 
the radius at which the binary is able to eject stars with 
greater than escape velocity. When there is a pre-existing 
core, the binary evolves somewhat differently than in 
these simulations; as we showed above, the orbital pe- 
riapsis progressively decreases while the apoapsis hardly 
changes. As a result, the orbital semi-major axis can still 
be large at the time that GW losses becomes signiflcant. 
Since most of the frictional force occurs near apoapsis, 
we do not expect significant stalling or core depletion to 
occur until late in the evol ution, perhaps not bef ore the 
two black holes merge (^e.g. lFukushige et al. Ill99"2[ ). Nev- 
ertheless, in what follows, we will explicitly note when in 
our integrations the mass of the sinking object becomes 
comparable to the mass in stars enclosed within its or- 
bital radius. 

We carried out calculations using the mass density pro- 
file of equation ([T4|) with {a = l]je = 1-8; 7 = 0.5; ro = 
600 pc;po = 35 Mqpc-^}, and M. = 3 x 10^ M© . The 
left panel of Figure [53| gives the orbital evolution of a test 
particle starting from an orbital radius of 100 pc, and ec- 
centricity Bin = (0.5,0.7,0.9). Dotted lines in the panel 
represent the radii at which the stellar mass enclosed in 
the orbit is 10'^ or 4 x 10^ Mq. For the two more ec- 
centric orbits (two left-most curves), it is possible that 
the binary enters the GW regime before violating these 
conditions. 

Although in our model the binary black hole mass is 
above the range (10^ — 10^ -^0) normally associated 
with space-based interferometers, we can nevertheless 
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Fig. 23. — Left panel: orbital evolution in the a, (1 — e) plane for a massive object in the M87 core starting from various eccentricities 
Ein = (0.5,0.7,0.9), and from an initial semi-major axis = 100 pc. Dynamical friction and gravitational-wave energy losses are both 
included. Dotted lines represent the radii at which the stellar mass enclosed in the orbit is 10^ M p, (lower curve) or 4 X 10^ Mq (upper 
curve). Red an d b lue lines are obtained respectively from the standard Chandrasekhar formula II18I I. which neglects fast-moving stars, and 
from equation l|31|l that assumes in addition a Maxwellian velocity distribution. Black curves are based on the more general equation l|40| l. 
Horizontal solid line gives the ISCO radius for a non-spinning hole (i.e., 6 gravitational radii). Central panel: time evolution of orbital 
semi-major axis (solid lines), apoapsis (upper-dashed curves) and periapsis (lower-dashed curves) in the three integrations with 6;^ = 0.7 
performed using: (i) the correct formula that includes the contribution from fast moving stars (black cu rve s), (ii) equation l|18| l in which 
only stars moving slower than the test mass contribute to the frictional drag (red curves); (iii) equation l|31|l which assumes a Maxewllian 
distribution of velocities (blue curves). Right panel: time evolution of apoapsis (solid lines) and orbital periapsis (dashed lines) for a 
4 X 10*^ Mq black hole. 



ask whether the eccentricity would remain large after 
the massive binary has entered into the GW regime. 
The Schwarzschild radius of a 3 x 10^ Mq SMBH is 
rsc ~ 1-4 X 10~''pc. When the orbital semi- major axis 
is 10 X rsc, we find that the binary eccentricity is still 
very large: e w (0.08, 0.6, 0.8) for ein = (0.5, 0.7, 0.9). 
When a = 5 X rsc the corresponding eccentricity is 
e « (0.03, 0.4, 0.7). 

The blue curves in Figure [23] were obtained by com- 
puting the frictional drag using Chandrasekhar's for- 
mula in its most common form, which assumes a locally 
Maxwellian distribution of velocities (equation [3T]). 
This approximation results in a very different orbital evo- 
lution characterized by smaller orbital eccentricities (for 
a given semi-major axis) and faster orbital decay when 
compared with the results obtained using the more cor- 
rect formula (PU)) . We note that - in spite of a higher rate 
of orbital decay - the smaller eccentricities achieved dur- 
ing the infall in this case result in a longer lifetime of the 
massive binary (central panel). The red curves in Fig- 
ure [23] were obtained using equation ()18|) , which allows 
for a non-Maxwellian velocity distribution but neglects 
the contribution to the frictional drag from stars moving 
faster than the sinking black hole. This approximation 
also results in a very different evolution when compared 
to the more correct treatment (black curve). Due to the 
smaller frictional drag, the standard treatment produces 
a slower decay of the orbital semi-major axis but a much 
faster evolution of the eccentricity which in turns results 
in a shorter life-time of the black hole binary. The right 
panel of Figure [23] shows the time evolution of orbital 
radius when M = 4 x 10^ Mq. In a shallow cusp near a 
SMBH, dynamical friction is very inefficient; this results 
in a very long sinking time. Starting from 100 pc, black 



holes with masses Af < 4 x 10^ Mq do not reach the 
center of the galaxy in a Hubble time unless their or- 
bit has a substantial initial eccentricity (ein < 0.7). We 
note however that such large eccentricities could be diffi- 
cult to retain at these radii due to orbital circularization 
that occurs outside the sphere of influence of the central 
SMBH. 

Cosmological simulations predict that a giant elliptical 
like M87 accreted about 4 Milky- Way-si zed galaxies over 
the last ~ 5 Gyr (jFakhouri et al. Il2010t ). The long sink- 
ing time scales found in Figure [53] suggest therefore that, 
at the present epoch, brightest cluster galaxies may still 
contain a few massive black holes or even satellite galax- 
ies (see below) moving through their extended cores. Al- 
though non-active secondary black holes could be very 
difficult to detect directly, such systems would be a pos- 
sible so urce of jet precession in the AGN of the central 
galaxy ([Romero et al. |[2000l) or they could induce a de- 
tectable displacement betwee n the galactic photo-cen ter 
and its nuclear point source ([Batcheldor et al. |[20lol ). 

In the computations presented above the infalling ob- 
ject was treated as a test particle of fixed mass. However, 
in a massive galaxy like M87 the central density is low 
enough that the infalling black hole may retain a signifi- 
cant fraction of stars from its host galaxy (because tidal 
forces are small). If stalling occurs, then one or more 
satellites may remain into the core of the central galaxy 
for a time significantly longer than a Hubble time. 

To address this possibility we integrated the equations 
of motion of a satellite galaxy in a fixed potential includ- 
ing the contribution of dyn amical friction and t he ef- 
fect of tidal truncation (e.g.. lAntonini et al. Il20lil ). The 
tidally truncated mass of the satellite galaxy (mr) is re- 
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Fig. 24. — Upper panel displays the orbital decay of satellite 
galaxies with different central velocity dispersions (or central black 
hole masses) into the core of M87. The evolution of the mass in 
stars of the infalling galaxies, as determined by the central galaxy 
tidal field, is given in the lower panel. 



lated to its limiting radius (ry) via 
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with a the one dimensional central velocity dispersion. 
The mass of the sa tellite SMBH is fixed by a through 
the M - cr relation (jGiiltiken et al. 1120091 ): 



M = 1.3 X 10*(cr/200 km s 
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The tidal radius can then be related to th e potential 
and density p of the central galaxy by fe.g.. lKiiig1ll96 
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Using for the central galaxy the mass distribution of 
equation p5|) . we find: 



dcf) Stt f r 

T = -T^Po'^o — 
dr \ro 



GM, 



(52) 



where po = 35 Mq/y>c^ and ro = 600 pc. This gives a 
limiting radius 



1 /9 n -1/2 

47r^ ^ 3GAf. 



(53) 



and a tidally-truncated mass from equation (|49p . Adopt- 
ing an initial distance of r = 600 pc and cr = 94 kms~^ 
(corresponding to Af = 4 x 10^ -AT©) we find rax = 
4.6 X 10'^ Mq and tt = 45 pc. 

Figure [M] plots the orbital evolution of satellites with 
initial orbital radius r — 600 pc and different values of 
the central velocity dispersion a = (50, 94, 200) km s^^ 
corresponding to M = (3 x 10^, 4 x 10'^, 10*^) Mq. In 
the core of a giant elliptical galaxy like M87, the time 
to reach the center for galaxies with a < 100 kms~^ is 
evidently longer than a Hubble time. 

First-ranked galaxies in clusters are often observed 
to contain multiple "nuclei," which may be identified 
wit h the tidally-tr uncated remains of inspiralling galax- 
ies (|Merritt 111981 . 

7. CONCLUSIONS 

In this paper we considered the orbital evolution of 
massive objects in nuclei with shallow density profiles 
around supermassive black holes (SMBHs). Our princi- 
ple results are summarized below. 

1 Orbital evolution can be very sensitive to the de- 
tails of the stellar distribution. In models with a 
flat central density profile, p ^ r^'', 7 w 0.5, the 
dynamical friction timescale is much longer than 
in models with a steep cusp due to the lack of low- 
velocity stars. The standard formula predicts that 
the inspiraling body will stall at a radius that is 
roughly 1/2 the core radius. 

2 Orbital eccentricity increases rapidly when the pe- 
riapsis falls inside the core. If the inspiralling body 
is initially at rbh with Cin ^ 0.5, its orbital eccen- 
tricity can become very large (> 0.9) by the time 
the orbit lies entirely inside the core. 

3 Using iV-body simulations, we found that the fric- 
tional force never falls precisely to zero. When the 
contribution of the fast-moving stars is included 
in the expression for the dynamical friction force, 
and (if appropriate) the changes induced by the 
massive body on the stellar distribution are taken 
into account, Chandrasekhar's theory reproduces 
the decay observed in the A^-body simulations very 
accurately. On the other hand, a straightforward 
application of Chandrasekhar's formula in its stan- 
dard form can give misleading results. 

4 If the mass of the inspiralling object is sufficiently 
large, it promotes the diffusion of stars into the 
phase-space region that was initially nearly empty, 
increasing the dynamical friction force. A low- 
density core is again regenerated as the object dis- 
places these stars. 

5 We derived an estimate of the Coulomb logarithm 
without any particular assumptions about the ve- 
locity distribution of field stars (e.g., that it follows 
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a Maxwellian distribution), and in the region out- 
side the core, where the standard dynamical fric- 
tion formula (jl8[) accurately represents the motion. 
We obtained In A w 6.5, consistent with previous 
work. 

6 We studied the location and evolution of the grav- 
itational wake that the inspiralling body induces 
in the stellar background. Outside the core, the 
peak in the overdensity lies close to the massive 
body at -20 < 6* < 0°, independent of M, and the 
amplitude of the overdensity increases with black 
hole mass. After the massive body enters the core, 
the density maximum decreases. This is consistent 
with the fact that the frictional drag is greatly re- 
duced inside the shallow cusp. 

7 In the absence of a steep central density cusp, the 
time for stellar-mass black holes to reach the cen- 
ter of the Milky Way from a starting radius of or- 
der 1 pc can easily exceed lOGyr. We computed 
the evolution of a population of stellar black holes 
as they segregate to the Galactic center, includ- 
ing the frictional force from the fast-moving stars. 
We found that, in models with parsec-scale cores, 
even after 10 Gyr, the density of black holes can 



remain substantially lower than that in stars at all 
radii. We conclude that it would be unjustified to 
assume that the massive remnants have yet reached 
a steady-state distribution at the Galactic center. 

8 Secondary black holes reach the gravitational radi- 
ation dominated regime on orbits that are typically 
very eccentric. However, we found that even ini- 
tially moderate eccentricities would result in non- 
negligible eccentricities at the moment the binary 
black hole enters the sensitivity window of planned 
space-based interferometers. This in turn would re- 
quire non-circular templates for gravitational wave 
data analysis. 

As a final remark, we recommend using equation p6p 
for the study of the inspiral of massive objects in galactic 
centers. 
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